function [R,F,Ab,LK] = resp_kernel(k, in)
% function [R,F,Ab,LK] = resp_kernel(keff, in)
% zeroth order response functions by semi analytical model

dat   = in.dat;
numg  = in.numg;
numm  = in.numm;
order = in.order;
if order > 0
    error(' order must be zero. ')
end

R    = zeros( in.numtypes, 4*numg*(order+1), 4*numg*(order+1) );
Rinc = zeros( numg*(order+1), numg*(order+1) );
Rfar = Rinc; Rlef=Rinc; Rrig=Rinc;
F    = zeros( in.numtypes, 4*numg*(order+1) );
Ab   = F;
tmpF = zeros( 1, numg*(order+1) );
tmpAb = tmpF;
LK = zeros( in.numtypes, 4,  numg*(order+1) );
FIS = zeros( 1, 2);
ABS = FIS;

% loop over all element types ---------------------------------------------
for el = 1:in.numtypes
    
    D1    = in.dat((el-1)*2+1, 1);
    D2    = in.dat((el-1)*2+2, 1);
    R1    = in.dat((el-1)*2+1, 2);
    A2    = in.dat((el-1)*2+2, 2);
    F1    = in.dat((el-1)*2+1, 3);
    F2    = in.dat((el-1)*2+2, 3);
    S12   = in.dat((el-1)*2+2, 5);
    Delta = in.xcm(2);
    
    % DO THE PHYSICS MODEL
    if (F2 > 0)
    
    %---------------------------------------------------------------------------
    % Modal Eigenvalues and related quantities
    
    % First eigenvalue (Maple-optimized; probably not good for the vector
    %   operations, but will translate to c)
    t1 = 1 ./ k;
    t4 = 1 ./ D1;
    t5 = (F1 .* t1 - R1) .* t4;
    t6 = 1 ./ D2;
    t7 = A2 .* t6;
    t9 = ((t5 + t7) .^ 2);
    t16 = sqrt((t9 + 4 * F2 .* t1 .* S12 .* t4 .* t6));
    t17 = t5 ./ 0.2e1 - t7 ./ 0.2e1;
    t18 = t16 ./ 0.2e1;
    lambda1 = t17 + t18;
    
    % Second eigenvalue (Maple-optimized)
    lambda2 = t17 - t18;
    
    % Modal matrix
    f11 = (F1./k-R1)./D1;
    f12 = F2./(D1.*k);
    f21 = S12./D2;
    f22 = -A2./D2;
    
    % FIT PARAMETERS FOR FULL ASSEMBLY with FISSION
    B=[    0.261513147949597
           0.014430493936852
          -1.963269877409992
          -1.686058019167120
           0.020446975313514
           0.211914109163937
           0.006730851120201
          -1.413949422121088
          -0.000148580985696
           0.612133597507877
           0.069427926629857
           0.562848227400144
           0.029221824260804
          -1.425806818665860
           0.000910188034673
           0.059643780247241
          21.962969406904719
          -2.368849865957258
          21.640984894797850
         -10.881107307101317];
     
      if (in.locfit==0)
          av1 = 0.0159;
          av2 = 0.0175;
          bv1 = 0.0015625.*Delta.^2 + 0.0125.*Delta.*D1 - 0.208333;
          bv2 = 0.0015625.*Delta.^2 + 0.0125.*Delta.*D2 - 0.208333;
      else
          av1 = B( 1)/Delta             + ...
              0.001*B( 2)             ;
          av2 = B( 3)/Delta             + ...
              0.02*B( 4)             ;
          bv1 = B( 5).*Delta + B( 6).*D1;
          bv2 = B( 7).*Delta + B( 8).*D1;   
      end
    
    v1 = sqrt(av1./bv1);
    v2 = sqrt(av2./bv2);

    %---------------------------------------------------------------------------
    % Effective modal eigenvalues and precomputed factors.
    kappa1 = lambda1 - v1.^2;% + beta(1)*D1 + beta(3)*D1*D2;
    kappa2 = lambda2 - v2.^2;% + beta(2)*D2 + beta(3)*D1*D2;
    
    K1 = (kappa1>=0).*sqrt(kappa1)   + (kappa1<0).*sqrt(-kappa1);
    C1 = (kappa1>=0).*cos(K1.*Delta) + (kappa1<0).*cosh(K1.*Delta);
    S1 = (kappa1>=0).*sin(K1.*Delta) + (kappa1<0).*sinh(K1.*Delta);
    K2 = (kappa1>=0).*sqrt(kappa2)   + (kappa1<0).*sqrt(-kappa2);
    C2 = (kappa1>=0).*cos(K2.*Delta) + (kappa1<0).*cosh(K2.*Delta);
    S2 = (kappa1>=0).*sin(K2.*Delta) + (kappa1<0).*sinh(K2.*Delta);
    
    %---------------------------------------------------------------------------
    % System of equations for a,b,c,d
    %   row 1
    m11 = 0.25.*(lambda1 - f22).*Delta;
    m12 = -0.5.*D1.*(lambda1-f22).*K1.*Delta;
    m13 = 0.25.*(lambda2 - f22).*Delta;
    m14 = -0.5.*D1.*(lambda2-f22).*K2.*Delta;
    %   row 2
    m21 = 0.25.*f21.*Delta;
    m22 = -0.5.*D2.*f21.*K1.*Delta;
    m23 = 0.25.*f21.*Delta;
    m24 = -0.5.*D2.*f21.*K2.*Delta;
    %   row 3
    m31 = (0.25.*(lambda1-f22)).*C1.*Delta - ...
        sign(kappa1).*(1/2).*D1.*(lambda1-f22).*S1.*K1.*Delta;
    m32 = (0.25.*(lambda1-f22)).*S1.*Delta+(0.5).*D1.*(lambda1-f22).*C1.*K1.*Delta;
    m33 = (0.25.*(lambda2-f22)).*C2.*Delta - ...
        sign(kappa1).*(0.5).*D1.*(lambda2-f22).*S2.*K2.*Delta;
    m34 = (0.25.*(lambda2-f22)).*S2.*Delta+(0.5).*D1.*(lambda2-f22).*C2.*K2.*Delta;
    %   row 4
    m41 = (0.25).*f21.*C1.*Delta - ...
        sign(kappa1).*(0.5).*D2.*f21.*S1.*K1.*Delta;
    m42 = (0.25).*f21.*S1.*Delta+(0.5).*D2.*f21.*C1.*K1.*Delta;
    m43 = (0.25).*f21.*C2.*Delta - ...
        sign(kappa1).*(0.5).*D2.*f21.*S2.*K2.*Delta;
    m44 = (0.25).*f21.*S2.*Delta+(0.5).*D2.*f21.*C2.*K2.*Delta;
    
    MM = [m11 m12 m13 m14; m21 m22 m23 m24; m31 m32 m33 m34; m41 m42 m43 m44];
    
    for g = 1:2
        x = 1; y = 0; x1 = [1;0;0;0];
        if g == 2
            x=0; y=1; x1 = [0;1;0;0];
        end
        
        [a, b, c, d] = invert(m11,m12,m13,m14, ...
                              m21,m22,m23,m24, ...
                              m31,m32,m33,m34, ...
                              m41,m42,m43,m44, x, y);
       y1 = MM \ x1; a = y1(1); b = y1(2); c = y1(3); d  = y1(4);
                          
                          
        %---------------------------------------------------------------------------
        % TOTAL GROUP FLUX RESPONSES
        
        phi1 = -Delta.*(b.*lambda1.*K2-b.*f22.*K2+d.*lambda2.*K1-d.*f22.*K1-...
            K2.*lambda1.*a.*S1-K2.*lambda1.*b.*C1+K2.*f22.*a.*S1+K2.*f22.*b.*C1-...
            K1.*lambda2.*c.*S2-K1.*lambda2.*d.*C2+K1.*f22.*c.*S2+K1.*f22.*d.*C2)./(K1.*K2);
        
        phi2 = f21.*Delta.*(-b.*K2-d.*K1+K2.*a.*S1+K2.*b.*C1+K1.*c.*S2+K1.*d.*C2)./(K1.*K2);
        
        %---------------------------------------------------------------------------
        % CURRENT RESPONSES -- OPTIMIZE THESE!
        
        % Group 1 out of Side 1
        R11g11 = (1./2).*D1.*(lambda2-f22).*K2.*Delta.*d + ...
            (1./4.*(lambda1-f22)).*Delta.*a + ...
            (1./4.*(lambda2-f22)).*Delta.*c + ...
            (1./2).*D1.*(lambda1-f22).*K1.*Delta.*b;
        % Group 2 out of Side 1
        R11g12 = (1./2).*D2.*f21.*K2.*Delta.*d+(1./4).*f21.*a.*Delta+...
            (1./4).*f21.*c.*Delta+(1./2).*D2.*f21.*b.*K1.*Delta;
        % Group 1 out of Side 2
        R12g11 = ((1./4.*(lambda1-f22)).*C1.*Delta-...
            (1./2).*D1.*(lambda1-f22).*S1.*K1.*Delta).*a+...
            ((1./4.*(lambda1-f22)).*S1.*Delta-...
            (1./2).*D1.*(lambda1-f22).*C1.*K1.*Delta).*b+...
            ((1./4.*(lambda2-f22)).*C2.*Delta-...
            (1./2).*D1.*(lambda2-f22).*S2.*K2.*Delta).*c+...
            ((1./4.*(lambda2-f22)).*S2.*Delta-...
            (1./2).*D1.*(lambda2-f22).*C2.*K2.*Delta).*d;
        % Group 2 out of Side 2
        R12g12 = ((1./4).*f21.*C1.*Delta-(1./2).*D2.*f21.*S1.*K1.*Delta).*a+...
            ((1./4).*f21.*S1.*Delta-(1./2).*D2.*f21.*C1.*K1.*Delta).*b+...
            ((1./4).*f21.*C2.*Delta-(1./2).*D2.*f21.*S2.*K2.*Delta).*c+...
            ((1./4).*f21.*S2.*Delta-(1./2).*D2.*f21.*C2.*K2.*Delta).*d;

        FIS(g)    = (phi1 .* F1         + phi2 .* F2).*Delta;
        ABS(g)    = (phi1 .* (R1 - S12) + phi2 .* A2).*Delta;
        
        % MORE FITTING
        if in.locfit == 1
        if g == 1
            R11g11 = R11g11 + B( 9)*Delta + B(10)*lambda1;
            R11g12 = R11g12 - 0.00482     + B(11)*lambda1;
            R12g11 = R12g11 + 0.006       + B(12)*lambda1;
            R12g12 = R12g12 + 0.00025     + B(13)*lambda1;
            ABS(1) = ABS(1) + B(17)*lambda1;
            FIS(1) = FIS(1) - 0.082 + B(19)*lambda1;
        else
            R11g11 = R11g11 - 0.043        + B(14)*lambda1;
            R11g12 = R11g12 + B(15)*Delta  + B(16)*lambda2;
            R12g11 = R12g11 + 0.003;
            % R12g12 = R12g12;         
            ABS(2) = ABS(2) + B(18)*lambda1;
            FIS(2) = FIS(2) + 0.038 + B(20)*lambda1;
        end
        end
        % fis/fis0 = a*p1*f1       + b*p2*f2
        % abs/abs0 = a*p1*(r1-s12) + b*p2*a2
        pp1 = (FIS(g)*A2-F2*ABS(g))/(F1*A2-R1*F2+S12*F2)/Delta;
        pp2 = (F1*ABS(g)-R1*FIS(g)+S12*FIS(g))/(F1*A2-R1*F2+S12*F2)/Delta;
        %k=1;
        % Balance correction     
        R13g11  = 0.5*(x - R11g11 - R12g11 + phi1*(F1/k-R1) + phi2*F2/k);
        R13g12  = 0.5*(y - R11g12 - R12g12 + phi1*S12       - phi2*A2  );
         
        % Balance correction     
        if (F2>0)
        R13g11  = 0.5*(x - R11g11 - R12g11 + pp1*(F1/k-R1) + pp2*F2/k );
        R13g12  = 0.5*(y - R11g12 - R12g12 + pp1*S12       - pp2*A2 );
        end
   
        if ( isnan(R13g11) )
            error('phi1 error')
        end
        if ( isnan(R13g12)  )
            error('phi1 error')
        end
        % Define the response matrices for this element
        Rinc(1,g) = R11g11;
        Rinc(2,g) = R11g12;
        Rfar(1,g) = R12g11;
        Rfar(2,g) = R12g12;
        Rlef(1,g) = R13g11;
        Rlef(2,g) = R13g12;
        if 1 == 0
            disp(['--- element: ',num2str(el)])
            fprintf(' R11g11:  %11.4f \n', R11g11)
            fprintf(' R11g12:  %11.4f \n', R11g12)
            fprintf(' R12g11:  %11.4f \n', R12g11)
            fprintf(' R12g12:  %11.4f \n', R12g12)
            fprintf(' R13g11:  %11.4f \n', R13g11)
            fprintf(' R13g12:  %11.4f \n', R13g12)
            fprintf('    FIS:  %11.4f \n', FIS)
            fprintf('    ABS:  %11.4f \n', ABS)
            fprintf('   phi1:  %11.4f \n', phi1)
            fprintf('   phi2:  %11.4f \n', phi2)
            disp('---')
            disp(' ')
        end
        
        LK(el,1,g) =  Delta*(R11g11 + R11g12 - 1); %incident
        LK(el,2,g) =  Delta*(R12g11 + R12g12    ); %far
        LK(el,3,g) =  Delta*(R13g11 + R13g12    ); %to right of incident
        LK(el,4,g) =  Delta*(R13g11 + R13g12    ); %to left of incident
        
    end
    % end groups
    
    % DO THE MODERATOR FITS
    else 
        
        [Rinc,Rfar,Rlef,FIS,ABS,LKb] = resp_kernel_mod(in);
        LK(el,:,:) = LKb;
    end
    
    
    Rrig      = Rlef;
    R(el,:,:) = [ Rinc  Rfar  Rrig  Rlef
                  Rfar  Rinc  Rlef  Rrig
                  Rlef  Rrig  Rinc  Rfar
                  Rrig  Rlef  Rfar  Rinc];
    F(el,:)   = [ FIS FIS FIS FIS]';
    Ab(el,:)  = [ ABS ABS ABS ABS]';
end
% end elements


end


function [a,b,c,d] = invert(...
    m11,m12,m13,m14, ...
    m21,m22,m23,m24, ...
    m31,m32,m33,m34, ...
    m41,m42,m43,m44, x, y)
% directly solve a 2x2 system

a = -(-m32.*m43.*x.*m24-m32.*m44.*m13.*y+m32.*m44.*m23.*x+m32.*m43.*m14.*y-m33.*m44.*m22.*x-m33.*m42.*m14.*y+m33.*m12.*m44.*y+m33.*m42.*x.*m24+...
    m34.*m42.*m13.*y-m34.*m42.*m23.*x+m34.*m22.*m43.*x-m34.*m12.*m43.*y)./(m23.*m41.*m32.*m14+m13.*m22.*m41.*m34-m13.*m41.*m32.*m24+m13.*m21.*m32.*m44+...
    m12.*m43.*m21.*m34+m12.*m41.*m33.*m24+m11.*m33.*m22.*m44-m21.*m32.*m43.*m14-m22.*m41.*m33.*m14-m21.*m33.*m12.*m44+m31.*m23.*m12.*m44-m23.*m11.*m32.*m44-...
    m23.*m12.*m41.*m34-m22.*m43.*m11.*m34+m11.*m32.*m43.*m24+m42.*m21.*m33.*m14-m42.*m13.*m21.*m34+m42.*m23.*m11.*m34-m31.*m42.*m23.*m14+m31.*m42.*m13.*m24+...
    m31.*m22.*m43.*m14-m31.*m13.*m22.*m44-m31.*m12.*m43.*m24-m42.*m11.*m33.*m24);
b = (-m31.*m43.*x.*m24-m31.*m44.*m13.*y+m31.*m44.*m23.*x+m31.*m43.*m14.*y-m44.*m21.*m33.*x+m34.*m21.*m43.*x+m41.*x.*m33.*m24+m44.*m11.*m33.*y-...
    m43.*m11.*m34.*y+m41.*m34.*m13.*y-m41.*m34.*m23.*x-m41.*m14.*m33.*y)./(m23.*m41.*m32.*m14+m13.*m22.*m41.*m34-m13.*m41.*m32.*m24+m13.*m21.*m32.*m44+...
    m12.*m43.*m21.*m34+m12.*m41.*m33.*m24+m11.*m33.*m22.*m44-m21.*m32.*m43.*m14-m22.*m41.*m33.*m14-m21.*m33.*m12.*m44+m31.*m23.*m12.*m44-m23.*m11.*m32.*m44-m23.*m12.*m41.*m34-m22.*m43.*m11.*m34+m11.*m32.*m43.*m24+m42.*m21.*m33.*m14-m42.*m13.*m21.*m34+m42.*m23.*m11.*m34-m31.*m42.*m23.*m14+m31.*m42.*m13.*m24+m31.*m22.*m43.*m14-m31.*m13.*m22.*m44-m31.*m12.*m43.*m24-m42.*m11.*m33.*m24);
c = -(m11.*m44.*y.*m32-m11.*m42.*m34.*y-m32.*m44.*m21.*x+m32.*m41.*x.*m24-m32.*m41.*m14.*y+m12.*m41.*m34.*y+m44.*m31.*m22.*x+m31.*m42.*m14.*y+m42.*m34.*m21.*x-m31.*m12.*m44.*y-m41.*x.*m22.*m34-m31.*m42.*x.*m24)./(m23.*m41.*m32.*m14+m13.*m22.*m41.*m34-m13.*m41.*m32.*m24+m13.*m21.*m32.*m44+m12.*m43.*m21.*m34+m12.*m41.*m33.*m24+m11.*m33.*m22.*m44-m21.*m32.*m43.*m14-m22.*m41.*m33.*m14-m21.*m33.*m12.*m44+m31.*m23.*m12.*m44-m23.*m11.*m32.*m44-m23.*m12.*m41.*m34-m22.*m43.*m11.*m34+m11.*m32.*m43.*m24+m42.*m21.*m33.*m14-m42.*m13.*m21.*m34+m42.*m23.*m11.*m34-m31.*m42.*m23.*m14+m31.*m42.*m13.*m24+m31.*m22.*m43.*m14-m31.*m13.*m22.*m44-m31.*m12.*m43.*m24-m42.*m11.*m33.*m24);
d = (m23.*m41.*m32.*x-m22.*m41.*m33.*x-m21.*m32.*m43.*x-m42.*m11.*m33.*y+m42.*m21.*m33.*x+m31.*m42.*m13.*y-m31.*m42.*m23.*x+m31.*m22.*m43.*x-m31.*m12.*m43.*y-m13.*m41.*m32.*y+m12.*m41.*m33.*y+m11.*m32.*m43.*y)./(m23.*m41.*m32.*m14+m13.*m22.*m41.*m34-m13.*m41.*m32.*m24+m13.*m21.*m32.*m44+m12.*m43.*m21.*m34+m12.*m41.*m33.*m24+m11.*m33.*m22.*m44-m21.*m32.*m43.*m14-m22.*m41.*m33.*m14-m21.*m33.*m12.*m44+m31.*m23.*m12.*m44-m23.*m11.*m32.*m44-m23.*m12.*m41.*m34-m22.*m43.*m11.*m34+m11.*m32.*m43.*m24+m42.*m21.*m33.*m14-m42.*m13.*m21.*m34+m42.*m23.*m11.*m34-m31.*m42.*m23.*m14+m31.*m42.*m13.*m24+m31.*m22.*m43.*m14-m31.*m13.*m22.*m44-m31.*m12.*m43.*m24-m42.*m11.*m33.*m24);

end
